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Abstract  We  compared  the  estimates  of  surface  drifter 
trajectories  from  1  to  7  days  in  the  equatorial  Atlantic 
over  an  18-month  period  with  five  eddying  ocean  gen¬ 
eral  circulation  model  (OGCM)  rcanalyscs  and  one 
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observational  product.  The  cumulative  distribution  of 
trajectory  error  was  estimated  using  over  7,000  days 
of  drifter  trajectories.  The  observational  product  had 
smaller  errors  than  any  of  the  individual  OGCM  rc¬ 
analyscs.  Three  strategies  for  improving  trajectory  es¬ 
timates  using  the  ensemble  of  five  operational  ocean 
analysis  and  forecasting  products  were  explored:  two 
methods  using  a  multi-model  ensemble  estimate  and 
also  spatial  low-pass  filtering.  The  results  were  insen¬ 
sitive  to  the  method  used  to  create  the  ensemble  esti¬ 
mates,  and  by  most  measures,  the  results  were  better 
than  the  observational  product.  Comparison  of  rela¬ 
tive  skill  of  the  various  OGCM  rcanalyscs  suggested 
promising  avenues  for  exploration  for  further  improve¬ 
ments:  forcing  with  higher  frequency  wind  stress  and 
quality  control  of  input  data.  One  of  the  lowest  hori¬ 
zontal  resolution  OGCMs,  with  1  /4°  longitude  horizon¬ 
tal  resolution,  made  the  best  trajectory  estimates.  The 
individual  OGCMs  were  dominated  by  errors  at  spatial 
scales  smaller  than  about  100  to  200  km,  i.e.,  less  than 
the  local  deformation  radius.  But  buried  in  those  errors 
were  valuable  signals  that  could  be  retrieved  by  com¬ 
bining  all  the  OGCM  velocity  fields  to  produce  a  multi- 
modcl  ensemble-based  estimate.  This  estimate  had  skill 
down  to  spatial  scales  about  75  km.  Results  from  this 
study  are  consistent  with  previous  work  showing  that 
ensemble-mean  forecast  skill  is  superior  to  individual 
forecasts. 

Keywords  Ocean  prediction  •  Surface  drifters  • 

Data  assimilation  •  Eddying  OGCM  • 

Model  intercomparison  •  Model-data  comparison  • 
Multi-model  ensemble  prediction 
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1  Introduction 

The  problem  of  predicting  or  estimating  trajectories  of 
objects  floating  near  the  sea  surface  arises  in  impor¬ 
tant  applications,  such  as  search  and  rescue  missions 
when  something  or  someone  is  lost  at  sea  at  a  known 
place  and  time,  or  when  debris  from  an  accident  is 
found  and  one  wishes  to  find  the  accident  location  by 
integrating  back  in  time  the  trajectories  of  the  debris. 
These  problems  require  knowledge  of  the  near  surface 
ocean  currents  and,  to  an  extent  depending  upon  the 
floating  object,  the  direct  effect  of  the  winds.  Here  we 
focus  on  the  first  requirement  of  accurate  near  surface 
currents  by  assessing  the  predictability  of  the  trajec¬ 
tories  of  Atlantic  Oceanographic  and  Meteorological 
Laboratory  (AOML)  surface  drifters  drogued  at  15  m 
depth.  We  focus  on  the  equatorial  Atlantic  since  this 
is  an  especially  challenging  area  due  to  the  presence 
of  strong  currents  having  a  large  space  and  time  vari¬ 
ability.  This  region  also  corresponds  to  the  area  where 
AF447  flight  crashed  (close  to  3°  N,  31°  W)  1st  June 
2009.  Several  research  laboratories  and  private  compa¬ 
nies  were  solicited  at  that  time  to  provide  estimates  of 
the  crash  position  given  the  location  of  the  recovered 
debris.  Using  ocean  currents  estimates,  the  backward 
trajectory  hindcast  of  these  floating  objects  would  be 
a  way  to  estimate  the  crash  location  (Drdvillon  et  al., 
this  issue).  This  was  a  good  motivation  to  study  surface 
drifter  trajectories  in  that  region. 

In  the  last  few  years,  operational  oceanography  cen¬ 
tres  have  developed  global  forecast,  nowcast  and  hind- 
cast  systems  providing  high  spatial  resolution  ocean 
currents  at  daily  or  higher  frequency.  These  data  are 
used  for  search  and  rescue  and  object  drift  applica¬ 
tions  (Davidson  ct  al.  2009)  as  well  as  marine  oil  pol¬ 
lution  prediction  (Hackctt  et  al.  2009).  Only  opera¬ 
tional  ocean  analysis  and  forecasting  systems  are  able 
to  predict  with  some  skill  the  ocean  surface  currents, 
especially  when  the  ocean  dynamics  is  dominated  by 
eddy  variability.  So  it  has  become  an  important  issue 
to  assess  not  only  surface  current  Eulcrian  statistics  but 
also  to  validate  the  Lagrangian  properties  of  the  ocean 
circulation  in  ocean  models.  The  capability  of  a  model 
to  forecast  a  drifter  trajectory  is  related  to  the  statistics 
of  the  distance  separating  two  particles  initially  close 
to  each  other.  One  expects  that  with  increasing  time, 
this  distance  will  increase  according  to  the  dynamics 
of  the  ocean  circulation.  The  study  of  the  separation 
distance  between  two  particles  was  first  addressed  by 
Richardson  (1926).  In  his  pioneering  study,  Richardson 
proposed  that  the  mean  squared  separation  distance 


<  d2  >  between  two  particles  is  proportional  to  /3 
(/  is  the  time).  This  behaviour  persists  until  d  ex¬ 
ceeds  the  size  of  the  largest  eddies  and  then  the 
separation  evolves  like  a  random  walk  process  and 

<  d2  >  becomes  proportional  to  time.  Thanks  to  qua- 
sigeostrophic  turbulence  and  surface  quasigeostrophic 
turbulence  theories,  it  has  been  possible  to  predict  the 
relative  dispersion  laws  as  a  function  of  the  slope  of  the 
kinetic  energy  wave  number  spectrum  (Sawford  2001; 
LaCasce  2008).  Numerous  studies  have  tried  to  ver¬ 
ify  these  theories  against  observations,  i.c.  balloon  (in 
the  atmosphere)  or  drifter  (in  the  ocean)  trajectories. 
LaCasce  (2010)  reviews  the  theoretical  probability  den¬ 
sity  functions  (PDFs)  used  to  characterize  pair  sepa¬ 
rations  and  tries  to  assess  their  validity  against  PDFs 
deduced  from  observations.  However,  pair  separation 
PDFs  are  difficult  to  validate  against  observations  as 
the  samples  considered  arc  often  small.  Results  also 
seem  to  be  sensitive  to  the  region  studied  as  noted 
by  LaCasce  (2008)  and  Lumpkin  and  Elipot  (2010). 
This  is  due  to  the  spatial  variability  of  the  velocity 
wave  number  spectrum.  Lagrangian  trajectories  are  in 
most  cases  isotropic.  However,  in  some  areas,  particle 
separation  may  be  affected  by  large-scale  mean  shear. 
In  the  case  of  a  zonal  front,  the  meridional  separation 
distance  will  be  different  from  the  zonal  one  (LaCasce 
2008).  Particles  can  also  be  affected  by  ocean  topogra¬ 
phy  (LaCasce  2008)  with  trajectories  tending  to  follow 
////  contours,  with  /  the  Coriolis  parameter  and  // 
the  ocean  depth.  Modelling  studies  have  investigated 
the  sensitivity  of  trajectory  forecasts  to  several  para¬ 
meters.  Ozgokmen  et  al.  (2000)  have  shown  that  given 
sufficient  surface  drifter  coverage,  RMS  forecast  errors 
of  less  than  50  km  can  be  obtained  for  periods  up  to 
3  months  in  the  gyre  interior  region  but  only  about 
1  week  in  western  boundary  current  and  midlatitudc  jet 
regions.  Their  study  also  shows  that  the  forecast  error 
is  very  sensitive  to  the  number  of  observations  avail¬ 
able  to  constrain  the  ocean  state  during  the  analysis. 
In  an  idealized  modelling  study,  Griffa  et  al.  (2004) 
revealed  that  the  spatial  smoothing  of  Eulerian  fields 
can  reduce  the  trajectory  forecast  skill.  Thanks  to  the 
availability  of  ocean  operational  forecasting  systems, 
studies  have  tried  to  estimate  float  trajectory  forecast 
skill  in  a  realistic  framework.  For  example,  Barron 
et  al.  (2007)  assessed  surface  drifter  trajectory  estimates 
in  many  regions  with  the  global  Navy  Coastal  Ocean 
Model  (NCOM)  model.  More  recently,  Huntley  et  al. 
(2011)  considered  trajectory  predictions  with  the  US 
Navy  East  Asian  Sea  (EAS16)  model  and  a  unique 
set  of  30  surface  drifters.  It  allowed  them  to  show 
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Table  1  Models  analysed  in  this  study 


Quantity 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

SURCOUF 

Eqns 

PE 

PE 

PE 

PE 

PE 

G,  Ek 

Mixing 

TKE 

TKE 

TKE 

MY2 

KPP 

Vert.  co. 

Z 

Z 

Z 

Z,<7 

p,  Z,  a 

n/a 

Nz 

50 

50 

50 

40 

(32) 

1 

dz  (m) 

2.25 

2.25 

2.25 

3.4 

2.4 

n/a 

dr (deg) 

1/4 

1/4 

1/12 

1/6 

1/12.5 

1/3 

Wind 

EC24 

EC24 

EC24 

NG3 

NG3 

ERA 

Output  ' 

Daily  mean 

Daily  mean 

Daily  mean 

6-hourly  snapshot 

Daily  snapshot 

6-hourly  mean 

TKE  is  Blanke  and  Delecluse's  (1993)  mixed  layer  scheme;  MY2  is  Mellorand  Yamada’s  (1974)  mixed  layer  scheme;  KPP  is  the  mixed 
layer  scheme  of  Large  et  al.  (1994) 

PE  primitive  equations,  <7  geostrophy,  Ek  Ekman,  Nz  total  number  of  levels  (layers),  dz  vertical  resolution  at  15  m,  dx  horizontal 
resolution  in  degree  longitude,  EC24  daily  ECMWF  stress,  ERA  six-hourly  ERA  Interim  stress,  NG3  three-hourly  NOGAPS 


that  predictive  skill  depended  more  upon  deployment 
location  than  time  of  deployment.  See  Huntley  et  al. 
(201 1 ,  and  references  therein)  for  a  summary  of  several 
other  related  studies. 

Several  questions  arise  in  using  operational  ocean 
forecasting  systems  to  estimate  drifter  trajectories  in 
the  real  ocean.  How  reliable  are  trajectory  forecasts 
made  in  a  realistic  framework?  Are  the  forecasts  very 
sensitive  to  the  choice  of  the  operational  surface  cur¬ 
rents  used?  Do  all  forecasts  have  the  same  skill  or  are 
there  differences  due  to  model  resolution,  data  assimi¬ 
lation  schemes  used  or  observations  assimilated?  Given 
an  ensemble  ocean  operational  forecasting  system, 
what  is  the  best  strategy  to  perform  the  most  accurate 


trajectory  forecast?  In  seeking  the  best  estimate,  should 
one  use  only  the  most  accurate  model?  Or  is  the  opti¬ 
mal  estimate  obtained  by  combining  the  estimates  from 
some  or  all  of  the  models?  How  should  one  combine  the 
model  estimates:  by  combining  their  velocity  fields  and 
performing  the  integration  with  the  compromise  veloc¬ 
ity,  or  by  performing  each  model  integration  separately 
and  combining  their  estimated  trajectory  end  points? 
We  addressed  these  questions  rigorously  by  assessing 
drifter  trajectory  estimate  success  using  products  sum¬ 
marized  in  Table  1  and  described  in  more  detail  in  the 
“Appendix”  and  in  Tables  2  and  3. 

Most  of  the  models  are  based  upon  primitive  equa¬ 
tion,  global  or  basin  scale,  ocean  general  circulation 


Table  2  Data  assimilated 


Data 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

Altimetry 

Jason- 1 

DM/NRT 

NRT 

NRT 

NRT 

NRT 

ENVISAT 

DM/NRT 

NRT 

NRT 

NRT 

NRT 

GFO 

DM/NRT 

NRT 

NRT 

NRT 

NRT 

SST 

Satellite 

NCEP 

NCEP 

NCEP 

NCODA 

NCODA 

RTG  0.5° 

RTG  0.5° 

RTG  0.5° 

SST 

SST 

AOML 

No 

No 

No 

Yes 

Yes 

drifters 

Subsurface 

T,  S 

ARGO 

DM  wQC 

NRT  wQC 

NRTwQC 

NRT  wQC 

NRTwQC 

Moored 

PIRATA, 

PIRATA, 

PIRATA, 

PIRATA, 

PIRATA, 

buoys 

wQC 

wQC 

wQC 

wQC 

wQC 

other 

XBT, 

XBT, 

XBT, 

XBT, 

XBT, 

CTD, 

CTD, 

CTD, 

CTD, 

CTD, 

MBATHY, 

MBATHY, 

MBATHY, 

MBATHY, 

MBATHY, 

etc  wQC 

etc  wQC 

etc  wQC 

etc  wQC 

etc  wQC 

NRT  near  real  lime,  DM  delayed  mode,  DM/NRT  before  Jan  23,  2008  DM  and  NRT  thereafter,  \\>QC  with  quality  control,  GFO 
Geosat  Follow-on,  NCODA  SST  see  Gentemann  et  al.  (2009) 
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Table  3  Data  assimilation  methods 

Method 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

IAU 

Yes,  7-davs 

No 

No 

Yes,  1-day 

Yes,  1-day 

Scheme 

Reduced 

Reduced 

Reduced 

NCODA 

NCODA 

order 

order 

order 

MVOI 

MVOI 

Kalman 

Kalman 

Kalman 

filter 

filter 

filter 

Velocity 

A 

B 

B 

C 

C 

initialisation 

A  statistical  by-product  of  the  analysis  (progressive  equatorial  cutoff  between  7°  N/7°  S),  B  physical  balance  operator,  C  NCODA 
MVOI,  geostrophically  balanced  increments  reduced  to  zero  at  2°  N/2°  S 


models  that  assimilate  observational  data.  The  excep¬ 
tion  is  SURCOUF,  which  uses  satellite  altimeter  data 
to  obtain  the  geoslrophic  velocity  and  wind  stress  to 
obtain  the  empirically  determined  Ekman  velocity.  We 
will  refer  to  all  these  velocity  products  as  “models'’ 
although  SURCOUF  is  more  of  an  observational  data 
product.  While  the  reliability  of  model  estimates  over 
various  integration  time  periods  is  an  important  test 
of  model  skill,  our  ranking  of  model  skill  was  only 
applicable  to  the  problem  at  hand;  for  example,  the 
models  with  higher-frequency  output  have  an  advan¬ 
tage  for  this  application  that  may  not  carry  over  to  other 
applications  that  do  not  require  high-frequency  output 
(for  a  comparison  of  metrics,  see,  e.g.  Metzger  et  al. 
2010). 

This  paper  is  organized  as  follows:  In  Section  2,  we 
delineate  the  study  area  and  in  Section  3,  we  describe 
the  methodology  used  to  estimate  surface  drifter  tra¬ 
jectories.  Then  in  Section  4,  we  present  the  trajectory 
estimate  skill  statistics  of  the  six  different  models.  In 
Section  5,  we  consider  how  to  improve  upon  the  esti¬ 
mates  of  the  single  OGCMs  by  combining  them  in  an 
optimal  way.  We  also  describe  the  results  of  applying  a 
spatial  filter  to  the  OGCM  velocity  fields,  which  both 
improved  the  estimates  and  allowed  us  to  interpret  the 
range  of  scales  for  which  the  models  had  useful  skill. 
We  conclude  with  a  summary  of  the  key  results  and 
their  implications  for  further  study  in  Section  6. 


2  Study  region 

The  comparison  of  various  surface  current  estimates  in 
the  tropical  Atlantic  was  first  motivated  by  the  need  to 
perform  drift  computations  in  the  context  of  the  search 
of  the  wreck  of  the  AF447  flight  from  Rio  to  Paris.  The 
airplane  disappeared  on  June  1st  2009  near  3°  N  and 
31°  W,  and  a  large  international  effort  was  organized  to 
try  to  find  the  wreckage.  The  contribution  of  oceanog¬ 
raphers  is  described  in  (Ollitrault  2010)  and  part  of  it  in 
Drdvillon  cl  al.  (2012).  This  region  is  very  challenging 


for  ocean  models  because  of  tropical  waves  and  es¬ 
pecially  for  SURCOUF  because  of  the  breakdown  in 
geostrophic  balance  at  the  equator.1  Here  we  might 
expect  the  dynamics  inherent  in  the  OGCMs  to  provide 
complementary  information  to  the  observational  data. 

The  choice  of  the  size  of  the  study  region  was  a 
subjective  compromise  between  choosing  a  larger  area 
that  included  more  surface  drifters  and  therefore  would 
provide  better  statistics  and  choosing  a  smaller  area 
that  was  more  homogeneous.  We  restricted  calculations 
to  drifters  within  the  main  study  region:  10°  S  to  10°  N, 
60°  W  to  0°  H,  outlined  in  Fig.  la.  We  studied  the 
18-month  period  from  July  2007  to  December  2008, 
the  period  when  the  models  in  Table  1  had  a  data 
assimilation  scheme  that  remained  unchanged. 

In  contrast  to  the  midlatitudcs,  this  region  has  mean 
flows  that  are  comparable  in  strength  to  that  of  the 
time-varying  flows.  Figure  la  shows  the  kinetic  energy 
(KE)  of  the  18-month  mean  flow  (MKE)  and  the  KE 
of  anomalies  relative  to  this  mean  (eddy  kinetic  energy, 
EKE)  for  the  SURCOUF  product.  Interestingly,  the 
mean  currents  reveal  more  sharply  defined  features 
associated  with  the  narrow,  zonally  aligned  current 
systems  of  the  equatorial  region.  For  example,  note  the 
MKE  minimum  from  about  35°  W  to  0°  at  about  2 
to  5°  N,  revealing  the  boundary  between  the  eastward 
flowing  north  equatorial  counter  current  (NECC)  to 
the  north  and  the  westward  flowing  south  equatorial 
current  (SEC)  to  the  south.  We  also  represented  in 
Fig.  lb  the  MKE  and  EKE  fields  for  the  ensemble 
multi-model  product  V5  described  in  Section  5.  The  V5 
product  is  a  combination  of  the  surface  velocity  fields 
from  the  five  OGCMs  used  in  this  study.  We  can  notice 
that  the  MKE  ofV5  (Fig.  lb)  is  much  smoother  and  less 
intense  than  in  SURCOUF  (Fig.  la).  The  EKE  field  of 


!More  properly,  the  meridional  velocity  becomes  ageostrophic 
while  the  zonal  geostrophic  velocity  becomes  balanced  with 
the  vertical  pressure  gradient  (Cushman-Roisin  and  Beckers 
2010),  which  cannot  be  diagnosed  from  sea  surface  height 
measurements. 
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(a)  log^MKE  [m2/s2])  :SURCOUF 


Fig.  1  a  Upper  panel  kinetic  energy  of  the  18-month  mean  flow 
from  the  SURCOUF  product.  Colour  bar  units  are  log,0(m2/s2). 
Low  er  panel  corresponding  kinetic  energy  of  the  current  anom¬ 
alies.  Large  rectangle  outlined  with  thin  straight  lines  indicates 
the  main  study  region.  Bold  polygon  outlines  the  equatorial 
subregion,  b  Like  in  panel  (a)  but  for  the  five-model  ensemble 
V5  product 


SURCOUF  and  V5  looks  also  quite  different.  V5  has  a 
stronger  EKE  along  the  equator.  This  may  be  related 
to  the  fact  that  OGCMs  have  daily  outputs  (NCOM 
outputs  are  even  six-hourly  snapshots)  whereas  SUR¬ 
COUF  surface  currents  comes  from  a  linear  interpola¬ 
tion  between  7-day  consecutive  altimetry  derived  cur¬ 
rents  maps.  We,  however,  can  sec  that  close  to  the 
Amazon  river  mouth,  SURCOUF  has  a  larger  EKE 
than  V5,  which  is  most  likely  due  to  tide  residuals  in  the 
SLA.  One  anticipates  that  it  will  be  important  for  mod¬ 
els  to  produce  both  accurate  time  variable  and  longer- 
term  mean  flows  for  accurate  trajectory  estimates. 


As  we  will  sec  later,  there  was  a  qualitative  difference 
between  the  sharpness  of  the  mean  flow  in  the  SUR¬ 
COUF  product  and  all  of  the  OGCM  mean  flows. 

In  contrast  to  the  sharply  defined  mean  flows,  the 
EKE  was  much  more  diffuse  (cf.  lower  panel  of  Fig.  la 
and  b).  T  his  is  due  to  the  equatorial  wave  guide  where 
linear  waves  (Rossby  waves  and  equatorial  Kelvin 
waves)  (Chelton  et  al.  2007;  Tulloch  ct  al.  2009)  travel 
along  the  equator.  A  relevant  feature  of  the  study 
region  is  the  occurrence  of  tropical  instability  waves 
(TIWs)  from  late  spring  through  to  January  (Weisbcrg 
and  Wcingartner  1988,  Fig.  6).  TIWs  are  surface- 
intensified  structures  of  600  to  1,200  km  wavelength 
with  phase  speeds  20  to  50  cm/s  and  current  speeds 
that  often  exceed  50  cm/s.  The  wave  periods  are  typ¬ 
ically  25  to  30  days  (Weisberg  et  al.  1987;  Cox  1980). 
TIWs  are  confined  to  the  equatorial  region,  with  EKE 
falling  off  very  rapidly  north  of  4°  N  (see  Jochum  ct  al. 
2004,  Fig.  7).  Despite  the  more  linear  dynamics  of  the 
tropics  and  wave  periods  of  about  a  month,  TIWs  can 
evolve  into  highly  nonlinear  tropical  instability  vortices 
(TIVs),  of  about  500  km  diameter  with  currents  reach¬ 
ing  over  1  m/s  (Foltz  et  al.  2004).  TIVs  have  been  ob¬ 
served  in  both  the  Pacific  (Kennan  and  Flament  2000) 
and  the  tropical  Atlantic  (Menkes  et  al.  2002;  Foltz 
ct  al.  2004).  OGCMs  resolving  the  equatorial  Rossby 
radius  arc  generally  able  to  simulate  the  observed  struc¬ 
tures  (Dutricux  et  al.  2008).  The  intcrannual  variability 
of  the  TIWs  is  difficult  to  capture  in  a  free-running  (no 
data  assimilation)  numerical  model  forced  with  realistic 
forcing  (von  Schuckmann  et  al.  2008). 

Another  notable  feature  of  the  study  region  is 
the  North  Brazil  Current  (NBC)  along  the  Brazilian 
coast.  This  current  rctroflccts  from  May  or  June  to 
December  and  feeds  the  North  Equatorial  Counter 
Current  (e.g.  Richardson  and  Rcverdin  1987).  The 
NBC  rclroflection  releases  NBC  rings  that  travel  along 
the  Brazilian  coast  towards  the  Lesser  Antilles  (Goni 
and  Johns  2001;  Barnier  et  al.  2001). 

The  large  TIV  intcrannual  variability  and  the  NBC 
rings  make  the  central  and  western  tropical  Atlantic  a 
region  where  currents  can  change  quickly.  Consider  for 
example  the  historical  time  scries  at  3°  N  and  28°  W 
from  the  SEQUAL  mooring  array  deployed  from  1983 
to  1985,  see  Fig.  2.  Around  day  201  of  1984,  the  current 
jumped  from  weakly  zonal  to  suddenly  very  strongly 
meridional  within  one  day,  see  closeup  in  the  lower 
panel  of  Fig.  2.  These  sudden  jumps  reveal  a  highly 
unstable  current  system  characterized  by  strong  inter¬ 
mittent  currents  reaching  over  100  cm/s.  For  that  par¬ 
ticular  event,  it  corresponds  to  a  very  intense  westward 
propagating  TI W  reaching  the  SEQUAL  surface  moor¬ 
ing  and  changing  abruptly  the  surface  current  intensity 
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Mooring  at  28W  3N,  Instrument  depth  10[m 


190  195  200  205  210 


day  of  the  year  1984 


Fig.  2  Northward  (black)  and  eastward  (red)  currents  from  the 
SEQUAL  surface  mooring  near  28°  W  3°  N.  Note  the  strong, 
chaotic,  highly  turbulent  currents  characteristic  of  this  barotrop- 
ically  unstable  region.  The  vertical  dashed  lines  bracket  a  20-day 
period  with  a  particular  strong  meridional  velocity  pulse.  Lower 
panel  shows  a  close-up  of  these  20  days  starling  on  day  190  of 
1984  (mid-July) 


(Weisberg  and  Weingartner  1988).  This  can  be  seen  on 
a  sequence  of  NOAA  1/4°  SST  maps  (Reynolds  et  al. 
2007,  not  shown). 

Others  have  noted  that  trajectory  estimate  skill 
strongly  depends  on  RMS  current  speed  (Ozgokmen 
et  al.  2001;  Barron  et  al.  2007).  Thus,  we  anticipated 
that  trajectory  estimate  skill  for  drifters  in  the  strong 
current  off  the  northwest  corner  of  Brazil  might  be  rela¬ 
tively  low,  while  skill  for  drifters  in  regions  with  weaker 
currents  (poleward  of  7°  N  and  4°  S)  might  be  relatively 
high.  Furthermore,  SURCOUF  velocities  were  com¬ 
puted  using  a  different  scheme  near  the  equator  (5°  S 
to  5°  N).  For  these  reasons,  we  defined  an  equatorial 
analysis  subregion  (see  Fig.  la).  Results  from  the  main 
study  region  and  the  equatorial  subregion  will  be  con¬ 
trasted  in  Section  4.2.  Results  in  Sections  4.1,  4.3  and  5 
were  for  the  main  study  region  only. 


3  Forward  trajectory  methodology 

We  used  the  observed  trajectories  of  satellite-tracked 
Surface  Velocity  Program  (SVP)  drifters  of  the  Global 
Drifter  Program  (http://www.aoml.noaa.gov/phod/dac/ 
gdp.html)  from  the  AOML  Drifter  Data  Assembly 
Center  (www.aoml.noaa.gov/phod/dac/dacdata.html). 
The  SVP  drifters  were  drogued  to  follow  the  currents 
near  15  m  depth,  with  the  ratio  of  the  drogue  drag  area 
to  remaining  drag  area  of  at  least  40.  The  wind-induced 
slippage  of  such  drifters  should  be  less  than  1  cm/s  in 
10  m/s  winds  according  to  Niiler  et  al.  (1995). 

Since  2005,  the  drifters  have  been  tracked  with  five 
or  six  satellites,  with  time  between  fixes  usually  be¬ 
tween  1  and  2  h  (Ellipot  and  Lumpkin  2008).  We 


used  the  quality-controlled  drifters  interpolated  to  six- 
hourly  positions  (Lumpkin  and  Pazos  2007). 

For  computational  convenience,  we  divided  the 
study  period  into  three  6-month  periods  from  July  1, 
2007  through  December  31,  2008.  All  the  drifter  tra¬ 
jectories  for  the  first  6-month  period  are  plotted  in 
Fig.  3.  For  each  drifter  in  the  study  region  and  for  each 
forecast  length  /i  e  {1,3,7}  days,  we  simulated  a  set  of 
trajectories  using  the  following  steps: 

1.  Start  at  its  earliest  location,  check  if  the  location  is 
known  n  days  later. 

2.  If  so,  estimate  this  later  location  using  second- 
order  Runge-Kutta  (Heun’s  method)  integration 
of  model  15  m  velocity  fields,  bilinearly  interpo¬ 
lated  in  space  and  time  and  using  the  model  out¬ 
put  time  step.  Trajectory  error  is  defined  as  the 
distance  (in  kilometres)  between  the  observed  and 
estimated  locations. 

3.  Check  that  each  model  simulated  trajectory  did 
not  run  aground.  If  this  occurs,  none  of  the  model 
forecast  is  taken  into  account  in  order  to  have  the 
same  number  of  forecasts  for  each  model. 

4.  Starting  from  the  end  time  of  the  last  forecast, 
repeat  the  above  procedure  (so  forecasted  trajec¬ 
tories  do  not  overlap  and  can  be  treated  as  indepen¬ 
dent).  If  the  start  and  end  location  of  a  drifter  are 
not  know,  reject  this  time  period  and  try  6  h  later 


Drifters  for  1-day  predictions  Jul-Dec  2007 


Fig.  3  Drifter  tracks  for  the  6-month  period  of  the  study  Jul 
1,  2007  to  Dec  31,  2007.  Red  symbols  are  the  drifters  outside 
the  equatorial  subregion.  Cyan  symbols  are  drifters  within  the 
equatorial  subregion,  which  is  also  outlined  by  the  bold  polygon 
(roughly  three  times  as  many  drifters  were  used  for  all  the  results 
of  this  study) 
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(i.c.  no  interpolation  was  done  of  the  six-hourly 
data). 

The  skill  of  each  model  is  then  estimated  by  comput¬ 
ing  the  median,  RMS  and  80th  percentile  of  the  trajec¬ 
tory  forecast  error  over  the  ensemble  of  forecasts  per¬ 
formed  (see  next  section).  We  indicated  the  confidence 
of  each  statistic  provided  by  a  bootstrap  method  (Efron 
1979).  We  also  computed  trajectories  backward  in  time. 
The  results  for  error  statistics  were  the  same  as  for 
the  forward  trajectories  to  within  statistical  uncertainty 
provided  by  the  bootstrap  method  so  we  have  only 
presented  forward  trajectories  herein. 


4  Single-model  trajectory  estimate  skill  results 

Here,  trajectory  estimate  skill  statistics  for  1-,  3-  and  7- 
day  trajectories  are  described  for  the  six  models  shown 
in  Table  1.  Results  are  shown  for  the  main  study  region 
(Section  4.1)  and  for  the  equatorial  subregion  (Sec¬ 
tion  4.2).  In  Section  4.3,  we  discuss  the  importance  of 
the  mean  flow  field. 

4.1  Main  study  region 

The  skill  of  each  model  is  summarized  in  the  cumula¬ 
tive  density  plot  of  trajectory  error  in  Fig.  4  and  the 
corresponding  statistics  in  Table  4.  The  80th  percentile 
of  trajectory  error  can  be  interpreted  as  the  radius  of 
the  circles  about  the  estimated  forward  trajectory  end 
point  that  contained  80%  of  the  actual  drifter  loca¬ 
tions.  We  could  have  shown  PDFs,  but  we  preferred 
to  present  the  results  with  cumulative  density  functions 
(CDFs)  because  it  is  easier  to  read  the  median  and 
80th  percentile  on  the  figure.  Furthermore,  a  CDF  is 
inherently  less  noisy  than  the  corresponding  PDF.  It  is 


Fig.  4  Cumulative  density  functions  of  trajectory  error  for  the 
main  analysis  region.  The  three  clusters  of  lines  from  left  to  right 
are  for  1-,  3-  and  7-day  estimates.  Solid  lines  are  for  the  six 
models,  labeled  by  their  first  letter  or  nvo,  e.g.  HYCOM,  shown 
in  green ,  is  labeled  “H”,  PSY2  is  labeled  “P2”.  The  comparison 
is  made  with  exactly  the  same  set  of  drifter  trajectories,  i.c.  if 
an  estimated  trajectory  ran  aground  or  otherwise  left  the  study 
region  for  one  model,  it  was  not  counted  for  any  of  the  models. 
Corresponding  values  tabulated  in  Table  4.  The  cyan  dashed  line , 
labeled  L5,  is  the  five-model  ensemble  of  OGCMs  (i.e.  without 
SURCOUF  )  combined  as  in  Eq.  1.  The  red  dashed  line ,  labeled 
V5,  is  the  five-model  ensemble  of  OGCMs  combined  as  in  Eq.  2 


worth  noting  that  the  forecast  error  PDFs  (not  shown) 
are  not  Gaussian.  They  resemble  the  separation  PDFs 
reviewed  by  LaCasce  (2010)  and  arc  strongly  asym¬ 
metrical.  The  PDF  peak  decreases  and  shifts  towards 
larger  separation  distances  with  increasing  trajectory 
time.  We  did  not  investigate  further  the  properties  of 
these  PDFs  as  it  is  beyond  the  scope  of  this  study. 

The  main  result  was  that  for  1-,  3-  and  7-day  tra¬ 
jectories,  SURCOUF  was  at  least  as  good  and  often 


Table  4  Single  model  forward-trajectory  error  and  95%  confidence  limits  for  the  main  study  region,  Jul  1, 2007  to  Dec  31, 2008 


Qntity 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

SURCF 

1-day 

median 

N  =  7,360 

16  ±0.3 

20  ±  0.4 

18  ±  0.4 

20  ±0.4 

20  ±0.4 

16  ±0.3 

RMS 

24  ±  0.5 

28  ±  0.6 

27  ±  0.6 

29  ±  0.6 

28  ±  0.5 

23  ±0.5 

80% 

28  ±  1 

32  ±  1 

31  ±  1 

34  ±  1 

34  ±  1 

27  ±  1 

3-day 

median 

N  =  2,365 

44  ±2 

50  ±2 

48  ±2 

53  ±2 

52  ±2 

38  ±  1 

RMS 

66  ±3 

75  ±3 

72  ±3 

82  ±3 

77  ±3 

60  ±2 

80% 

77  ±3 

86  ±3 

82  ±3 

93  ±4 

92  ±  4 

70  ±  3 

7-day 

median 

N  =  933 

93  ±6 

104  ±6 

98  ±6 

114  ±6 

111  ±8 

80  ±  4 

RMS 

147  ±  10 

160  ±  11 

149  ±9 

175  ±  10 

166  ±  10 

128  ±  7 

80% 

164  ±  10 

183  ±  12 

167  ±  10 

205  ±  12 

196  ±  12 

145  ±  10 

80%  means  80th  percentile.  Results  shown  graphically  in  Fig.  4 
N  the  number  of  estimates 
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much  better  than  the  much  more  expensive  OGCMs 
that  all  assimilated  altimeter  data.  Of  the  OGCMs, 
Global  Ocean  Reanalysis  and  Simulations  (GLORYS) 
made  the  best  estimates,  and  NCOM  made  the  worst 
estimates.  Prototype  Syst£me  2  (PSY2)  and  Prototype 
Syst£me  3  (PSY3)  had  similar  skill  though  PSY2  was 
often  slightly  better.  It  seems  that  the  quality  of  obser¬ 
vations  and  the  method  used  to  assimilate  them  may 
be  more  important,  since  the  GLORYS  trajectories 
were  the  best  among  all  OGCMs  despite  the  fact  that 
the  GLORYS  (and  PSY3)  grid  also  had  the  coarsest 
resolution.  This  highlights  the  importance  of  the  data 
quality  as  GLORYS  is  a  reanalysis  using  reprocessed 
and  quality  checked  datasets  whereas  PSY2  and  PSY3 
operational  outputs  come  from  the  real-time  operation 
of  the  systems.  The  mode  of  data  assimilation  also  has 
an  impact  as  a  smooth  initialisation  with  incremental 
analysis  update  is  used  in  GLORYS  instead  of  sequen¬ 
tial  correction  in  (higher  resolution)  PSY2. 

The  forecast  skill  of  a  simple  random  walk  model 
was  also  tested  to  see  if  the  trajectory  forecast  in  the 
models  used  would  be  similar  to  a  diffusive  process.  To 
do  this,  we  used  the  approach  described  in  LaCasce 
(2008)  (see  their  Eq.  23).  The  random  walk  model 
parameters  are  the  local  time  average  and  root  mean 
square  of  the  velocity.  In  our  case,  we  used  the  velocity 
fields  from  GLORYS.  The  results  (not  shown)  indicate 
that  the  forecast  skill  of  the  random  walk  model  is 
significantly  worse  than  GLORYS.  This  shows  that  at 
the  scales  considered  (1-,  3-  and  7-day  forecasts),  the 
model  currents  are  not  a  simple  noise. 

4.2  Equatorial  subregion 


Fig.  5  As  in  Fig.  4  but  for  the  equatorial  subregion.  Correspond¬ 
ing  single-model  values  tabulated  in  Table  5 


corresponding  statistics  summarized  in  Table  5.  Results 
for  this  subregion  arc  very  similar  to  those  for  the  main 
analysis  region.  In  fact  the  trajectory  errors  are  often 
slightly  smaller  for  the  equatorial  subregion,  though 
the  differences  are  not  significant  at  the  95%  level. 
With  only  about  1/2  as  many  estimates  (because  of  the 
smaller  region),  the  statistical  uncertainty  was  larger. 
In  general,  the  model  trajectory  errors  were  similar. 
Note  especially  that  SURCOUF  consistently  made  the 
best  estimates  in  the  subregion  as  well,  with  the  minor 
exception  of  the  extreme  tails  of  the  distributions,  cf. 
Fig.  5. 

The  similar,  or  perhaps  slightly  better,  trajectory  es¬ 
timates  in  the  equatorial  subregion  are  both  surprising 
and  informative.  Recall  the  elevated  EKE  associated 


Cumulative  density  plots  of  trajectory  error  for  the  with  TIWs  fell  off  quickly  northward  of  4°  N  and  south- 
equatorial  subregion  are  shown  in  Fig.  5,  with  the  ward  of  1°  S  (Jochum  et  al.  2004,  Fig.  7),  and  TIVs  form 


Table  5  Like  Table  4  but  for  the  equatorial  subregion  outlined  in  Fig.  la 


Quantity 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

SURCF 

1-day 

median 

N  =  3,864 

17  dbO 

21  ±0 

19  db  0 

20  db  0 

22  db  1 

16  dbO 

RMS 

23  i  1 

27  ±  1 

26  ±  1 

27  ±  1 

28  db  1 

22  db  0 

80% 

28  ±  1 

32  db  1 

31  ±  1 

34  db  1 

35  ±  1 

27  db  1 

3-day 

median 

N  -  1 ,226 

48  ±2 

55  ±  3 

52  db  2 

57  ±3 

58  db  2 

41  ±2 

RMS 

64  ±3 

72  ±3 

70  ±  3 

76  ±3 

79  ±3 

58  db  2 

80% 

78  ±3 

86  db  3 

84  db  4 

94  db  4 

97  ±5 

72  ±4 

7-day 

median 

N  —  477 

99  ±9 

I10±9 

104  ±7 

121  ±9 

121  ±9 

85  ±7 

RMS 

140  db  13 

151  db  14 

1 50  db  1 3 

163  db  9 

1 70  db  1 0 

124  db  8 

80% 

167  ±  11 

175  db  10 

176±  12 

205  ±  15 

209  db  19 

144±  12 

The  statistical  uncertainty  is  greater  because  of  the  smaller  numbers  of  estimates 
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along  the  high  lateral  shear  zone  between  the  NECC 
and  SEC  and  are  mostly  found  between  the  equator 
and  about  5°  N  (Foltz  et  al.  2004,  Fig.  5a).  Thus,  these 
challenging  features  are  mostly  within  our  equatorial 
subregion,  and  thus  cannot  be  the  sole  reason  limiting 
the  performance  of  the  OGCM  trajectory  estimates  nor 
the  sole  reason  for  poor  OGCM  performance  relative 
to  the  SURCOUF  data  product.  This  uniform  perfor¬ 
mance  of  SURCOUF  was  surprising  since  gcostrophic 
balance  “breaks  down”  in  the  equatorial  region,  and 
close  to  the  equator,  only  the  zonal  velocity  component 
of  the  velocity  is  deduced  from  SLA. 

4.3  18-month  mean  flow 

Given  the  relatively  strong  mean  flows  in  the  tropical 
Atlantic,  a  natural  question  is  whether  the  18-month 
mean  flow  of  SURCOUF  was  perhaps  superior  to 
that  of  the  OGCMs.  In  Figs.  6  and  7,  we  plot  the 
mean  zonal  and  meridional  velocity  fields,  respectively, 
for  several  of  the  models.  All  three  Mcrcator-Occan 
OGCM  (i.c.  GLORYS,  PSY3,  PSY2)  mean  fields  were 
so  similar  that  they  were  difficult  to  distinguish  visually 
and  thus  only  GLORYS  is  presented.  While  there  are 
slight  differences  between  GLORYS  and  Hybrid  Co¬ 
ordinate  Ocean  Model  (HYCOM)  in  both  fields,  the 
most  unique  field  is  the  zonal  velocity  of  SURCOUF, 
with  sharply  defined  structures  such  as  the  SEC.  To  test 
whether  both  the  mean  and  the  eddy  fields  of  SUR¬ 
COUF  had  more  skill  than  the  corresponding  OGCM 
fields,  we  performed  in  the  main  region  several  ex¬ 
periments  in  which  we  replaced  the  SURCOUF  mean 


Fig.  6  Eighteen-month  mean  zonal  velocity  field  for  a  SUR¬ 
COUF,  b  GLORYS  and  c  HYCOM 


(a) 


fields  with  those  from  other  models.  In  all  cases,  the 
trajectory  errors  increased  although  not  always  enough 
to  be  statistically  significant  at  the  95%  level.  This  was 
true  even  for  the  V5  ensemble  (defined  in  Section  5). 
This  confirmed  that  the  fine  structures  in  Fig.  6a  were 
an  improvement  over  the  more  diffuse  mean  zonal 
velocities  of  the  OGCMs. 

One  can  also  ask  whether  the  mean  field  of  SUR¬ 
COUF  was  the  most  important  factor?  That  is,  perhaps 
the  eddy  field  of  SURCOUF  was  actually  inferior  to  the 
OGCMs,  but  the  deficiency  was  compensated  for  by  the 
superior  mean  field?  We  confirmed  that  this  was  not  the 
case.  We  replaced  the  GLORYS  mean  field  with  that 
of  SURCOUF  and  found  that,  while  this  improved  the 
trajectory  estimates  over  that  of  GLORYS  alone,  the 
estimates  were  still  inferior  to  SURCOUF  alone.  That 
is,  the  eddies  of  SURCOUF  were  also  superior  to  the 
eddies  of  GLORYS.  However,  as  we  will  see  in  the  next 
section,  the  ensemble  average  of  all  the  OGCMs,  V5, 
made  better  trajectory  estimates  than  SURCOUF.  This 
occurred  despite  the  fact  that  its  mean  field  was  inferior 
to  that  of  SURCOUF.  Thus,  the  V5  ensemble  seems 
to  have  superior  eddies  than  SURCOUF,  although  this 
could  not  be  explicitly  verified. 


5  Improving  model  trajectory  estimates 

We  might  expect  that  somehow  combining  models 
might  improve  the  trajectory  estimates.  The  benefits 
from  ensemble  forecasts  on  predictability  has  been  no¬ 
ticed  first  by  Leith  (1974),  more  than  three  decades  ago. 
The  idea  behind  ensemble  forecasting  is  that  an  ensem¬ 
ble  of  analyses  or  forecasts  will  better  sample  the  true 
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state  than  a  single  member  (deterministic  forecast). 
By  averaging  the  ensemble  members,  one  gets  an  im¬ 
proved  true-state  estimation  than  with  a  single  member. 
A  comprehensive  description  of  the  ensemble  forecast¬ 
ing  concept  is  given  by  Murphy  (1988).  More  recently, 
the  “super  ensemble”  approach  of  optimally  combining 
several  different  model  estimates  has  been  successfully 
used  in  ocean  forecasting  (Rixen  and  Ferreira-Coelho 
2007;  Rixen  et  at.  2008,  2009).  However,  for  the  ap¬ 
plication  explored  herein,  the  best  method  to  combine 
models  is  not  a  priori  obvious.  On  the  one  hand,  one 
could  first  use  each  model  to  make  an  estimate  and 
then  combine  those  estimates.  The  advantage  with  such 
as  method  is  that  each  model  is  internally  dynamically 
consistent.  On  the  other  hand,  one  could  combine  the 
model  velocity  fields  and  make  estimates  with  these 
compromise  velocity  Fields.  The  apparent  advantage 
of  such  a  method  is  that  trajectory  errors  accumulate 
during  the  trajectory  integration  so  it  might  be  best  to 
reduce  errors  in  the  velocity  fields  before  using  them 
for  integration.  Furthermore,  the  compromise  velocity 
fields  should  be  smoother,  reducing  the  noise  we  expect 
to  dominate  at  small  scales.  Below  we  describe  the 
results  from  employing  these  two  strategies.  We  also 
address  the  possibility  of  improving  the  single-model 
estimates  by  applying  a  low-pass  spatial  filter  to  the 
velocity  fields.  This  analysis  allowed  us  to  identify  the 
range  of  length  scales  with  useful  skill  for  both  the 
individual  OGCMs  and  the  five-model  compromise. 

5.1  Four-  vs.  five-model  ensembles 

We  computed  compromise  forecasts  from  an  ensem¬ 
ble  of  four  or  five  models.  Consider  the  “model 


compromiseforecast”,  as  the  weighted  average  of  the 
estimated  locations: 


.yens  —  I  VJffl  (0 

where  xm  and  ym  arc  the  estimated  longitude  and  lat¬ 
itude  of  the  forward  trajectory  integrations  for  model 
m.  For  weights  we  used,  wm  —  i  where 

Vm  was  the  mean  square  error  for  model  m  for  the  cor¬ 
responding  forecast  length.  For  the  ensemble  of  all  five 
OGCM  models,  the  sum  is  over  m  e  { 1 , 2, . . .  5  =  A/}. 

Statistics  for  the  trajectory  errors  corresponding  to 
the  five-model  ensemble  are  presented  in  the  seventh 
column  of  Table  6,  labeled  “ALL”.  We  also  considered 
removing  one  model  from  the  ensemble.  The  results  for 
these  four-model  ensembles  are  also  presented  in  Ta¬ 
ble  6  in  columns  2  through  6.  The  column  label  indicates 
the  model  withheld  from  the  ensemble.  Note  that  none 
of  the  columns  has  trajectory  error  consistently  lower 
than  the  seventh  column.  This  indicates  that  none  of  the 
models  consistently  degraded  the  five-model  ensemble 
estimate.  All  models  “contributed”  in  the  sense  that 
including  them  in  the  five-model  ensemble  made  for  a 
better  estimate  or  at  least  not  a  worse  estimate. 

We  expect  the  improvement  from  the  multi-model 
ensemble  to  arise  from  the  cancellation  of  random 
model  trajectory  errors.  Ideally  all  the  models  would 
provide  independent  estimates,  so  their  errors  would  be 
unrelated  and  tend  on  average  to  cancel.  The  cumula¬ 
tive  density  function  for  the  five-model  OGCM  ensem¬ 
ble,  composed  of  GLORYS,  PSY3,  PSY2,  NCOM  and 
HYCOM  combined  as  in  Eq.  1,  was  plotted  in  Fig.  4. 
This  compromise  estimate  composed  of  only  OGCM 


Tabic  6  OGCM  model  ensemble  forward  trajectory  error  and  95%  confidence  limits  for  the  main  study  region  for  time  period  July  1, 
2007  to  Dee  31, 2008 


Qntity 

GLORYS 

PSY3 

PSY2 

NCOM 

HYCOM 

ALL 

1-day 

median 

N  =  7,360 

16  ±0.3 

15  ±0.3 

15  ±0.3 

15  ±0.3 

15  ±0.3 

15  ±0.3 

RMS 

23  ±  0.4 

22  ±  0.4 

22  ±0.4 

21  ±0.4 

22  ±  0.4 

21  ±0.4 

80% 

27  ±0.5 

26  ±0.5 

26  ±  0.5 

26  ±  0.4 

26  ±  0.5 

26  ±  0.4 

3-day 

median 

N  =  2,365 

41  ±  1 

39  ±2 

40  ±  1 

39  ±  1 

40  ±  1 

39  ±  1 

RMS 

60  ±2 

58  ±2 

58  ±2 

57  ±2 

59  ±2 

57  ±2 

80% 

71  ±3 

68  ±2 

69  ±2 

68  ±2 

69  ±2 

68  ±2 

7-day 

median 

N  =  933 

88  ±4 

84  ±5 

84  ±  4 

84  ±3 

85  ±5 

82  ±5 

RMS 

127  ±7 

125  ±8 

124  ±  6 

122  ±7 

126  ±7 

122  ±  7 

80% 

149  ±8 

145  ±9 

149  ±7 

142  ±  9 

149  ±8 

143  ±7 

The  seventh  column  “ALL”  is  for  the  weighted  average  of  all  five  OGCM  models.  The  remaining  columns,  second  through  sixth,  are  for 
the  average  of  four  models.  So  the  second  column  is  for  all  but  the  GLORYS  model.  The  higher  trajectory  error  for  the  second  column 
indicates  that  removing  GLORYS  from  the  ensemble  worsens  the  ensemble  estimate.  None  of  the  models  were  found  to  consistently 
worsen  the  estimates 
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results  was  clearly  as  good  as  or  sometimes  better  than 
the  data  product  SIJRCOUF  and  all  the  individual 
model  estimates.  This  result  further  supports  the  idea 
of  partial  cancellation  between  model  trajectory  errors. 

5.2  Combining  model  velocity  fields 

Another  strategy  for  combining  model  estimates  is  to 
average  the  velocity  fields  from  an  ensemble  of  models 
and  make  a  single  integration  with  the  ensemble  veloc¬ 
ity  field.  Consider  the  “model  compromise  velocity”,  as 
the  weighted  average  of  the  velocity  from  the  individual 
models: 


Uens  —  5Zw=l  Vm  (2) 

where  um  and  vm  are  the  zonal  and  meridional  ve¬ 
locities  for  the  15-m  currents  for  model  m,  linearly 
interpolated  onto  a  common  grid.  For  weights  we  used, 
wm  =  F"1/ V~\  where  Vm  was  the  mean  square 
error  for  model  m  for  the  7-day  trajectory  forecasts. 
For  the  common  grid,  we  used  that  of  the  GLORYS 
model.  Note  that  this  effectively  coarsened  the  spatial 
resolution  of  all  models  except  PSY3  and  coarsened  the 
temporal  resolution  of  NCOM.  We  expect  the  spatial 
coarsening  had  negligible  effect  based  upon  the  results 
of  our  smoothing  experiments  reported  in  the  next 
subsection.  We  confirmed  that  the  temporal  smooth¬ 
ing  had  negligible  effect  by  repeating  the  experiment 
using  the  NCOM  temporal  resolution  (six  hourly)  and 
found  agreement  well  within  statistical  uncertainty. 
Hereinafter  we  will  refer  to  these  fields  as  “V5”. 

The  empirical  cumulative  density  function  of  trajec¬ 
tory  error  corresponding  to  V5  ensemble  was  plotted  in 
Fig.  4  with  a  dashed  red  line.  Interestingly,  both  strate¬ 
gies  for  combining  the  model  estimates  gave  similar 
results  (both  red  and  cyan  dashed  lines  are  similar). 
And  regardless  of  the  strategy  to  combine  the  model 
estimates,  the  five-model  OGCM  ensembles  (super  en¬ 
sembles)  both  worked  better  than  all  the  individual 
OGCMs  and  had  similar  skill  and  sometimes  better  skill 
than  the  data  product  SIJRCOUF. 

5.3  Smoothing  velocity  fields 

One  might  speculate  that  combining  the  velocity  fields 
as  in  Eq.  2  led  to  smoother  velocity  fields.  Pessimisti¬ 
cally  one  might  conjecture  that  the  improved  estimates 
occurred  only  because  of  the  smoothing  (Wunsch  2010, 
personal  communication),  a  position  we  will  call  the 
null  hypothesis.  More  generally,  it  seems  plausible  that 
the  errors  in  the  OGCM  velocity  fields  were  a  strong 


function  of  spatial  scale  with  larger  errors  at  smaller 
scales.  This  raises  the  question  as  to  what  length  scales 
have  useful  skill  and  what  scales  are  too  small  to  be  cap¬ 
tured  by  the  OGCMs.  To  address  this  question,  we  have 
applied  a  low-pass  filter  to  the  OGCM  velocity  fields 
prior  to  Runge-Kutta  trajectory  integration.  The  low- 
pass  filter  was  accomplished  by  convolving  the  original 
OGCM  velocity  fields  at  each  time  step  and  each  grid 
point  (jt0,  yo)  with  a  Gaussian  filter, 

G(x,  y;  x0,  y0)  =  A  exp{-[(*  -  x0)2  +  (y  -  yo)2V(2R2)}, 

with  Gaussian  radius  /?,  over  a  square  region  centred 
at  (*0,yo)  and  of  width  2/?,  with  A  chosen  so  that 
the  sum  of  weights  was  unity.  For  each  OGCM,  many 
experiments  were  performed  with  R  ranging  from  a 
few  tens  of  kilometres  to  hundreds  of  kilometres,  each 
experiment  using  only  one  application  of  a  given  filter 
to  the  original  velocity  fields. 

There  appeared  to  be  an  optimal  smoothing  radius 
in  the  range  100  <  R  <  200  km  that  produced  the  best 
estimates.  This  is  shown  with  the  80th  percentile  of  the 
7-day  trajectory  error  plotted  in  Fig.  8  versus  smoothing 
radius  R.  The  corresponding  plot  for  the  RMS  error 
is  shown  in  Fig.  9.  The  y-intercept,  corresponding  to 
R  =  0,  corresponds  to  the  result  with  no  smoothing  at 
all,  for  which  values  can  also  be  read  from  the  final  two 
rows  of  Table  4,  though  agreement  is  not  exact  since 


Effect  of  spatial  smoothing  on  prediction  error 


Fig.  8  Eightieth  percentile  of  7-day  trajectory  error  (kilometres) 
vs.  smoothing  radius  (kilometres)  for  each  OGCM.  Analysis 
performed  in  main  study  region.  Line  labelling  described  in  Fig.  4 
caption.  Error  bars  represent  the  formal  error  95%  confidence 
limits.  Discrepancies  between  the  K-intercept  values  and  the  final 
row  of  Table  4  arise  because  different  trajectories  were  used. 
Here  near  coastal  trajectories  were  eliminated  when  the  velocity 
fields  were  smoothed 
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Effect  of  spatial  smoothing  on  prediction  error 


Fig.  9  As  in  Fig.  8  but  for  the  RMS  error.  Discrepancies  between 
the  /-intercept  values  and  the  penultimate  row  of  Table  4  arise 
because  different  trajectories  were  used.  Here  near  coastal  tra¬ 
jectories  were  eliminated  when  the  velocity  fields  were  smoothed 


trajectories  near  the  coasts  were  necessarily  omitted 
from  Figs.  8  and  9,  due  to  the  smoothing.  In  fact  the 
discrepancies  underline  the  importance  of  always  mak¬ 
ing  comparisons  with  the  same  set  of  drifter  trajectories 
because  of  the  long  tails  in  the  PDFs.  Error  bars  rep¬ 
resent  the  95%  confidence  interval  obtained  with  bool 
strapping  with  1,000  subsamplcs.  Because  of  the  large 
statistical  uncertainly,  it  is  necessary  to  consult  both 
Figs.  8  and  9  to  confirm  the  trends.  The  local  minimum 
in  trajectory  error  was  qualitatively  reproduced  in  both 
statistics,  with  the  exception  of  the  80th  percentile  for 
PSY2. 

The  interesting  exception  to  the  optimal  smoothing 
radius  was  the  model  constructed  from  all  five  OGCMs 
as  per  Eq.  2,  labeled  V5  in  Figs.  8  and  9,  which  showed 
little  or  no  improvement  for  small-scale  smoothing, 
and  errors  worsened  noticeably  for  R>  75  km.  The 
errors  continued  to  grow  monotonically  and  rapidly 
through  to  the  largest  filter  radius  considered,  R  % 
450  km.  That  is,  the  optimal  smoothing  radius  was 
indistinguishable  from  zero  (no  smoothing)  and  was  at 
most  about  75  km.  Unfortunately  it  was  difficult  to  be 
more  precise  because  the  95%  confidence  limits  were 
fairly  large.  Accounting  for  the  statistical  uncertainty, 
one  can  say  the  V5  results  in  Figs.  8  and  9  show  with 
95%  confidence  that  the  trajectories  with  about  150  km 
or  larger  smoothing  were  worse  than  those  with  no 
smoothing.  This  implies  that  there  is  at  least  some 
skill  in  length  scales  less  than  150  km.  However,  the 
smoothness  of  the  curves  and  the  consistency  of  the 
80lh  percentile  of  the  7-day  trajectory  error  and  RMS 


error  suggests  these  arc  conservative  confidence  limits. 
The  smaller  optimal  smoothing  radius,  between  0  and 
75  km,  for  V5  compared  to  100  to  200  km  for  any  of  the 
individual  OGCMs  suggests  that  V5  had  skill  at  smaller 
scales.  Our  interpretation  is  that  buried  in  the  noise 
at  small  scales,  the  individual  OGCMs  had  some  skill 
down  to  about  75  km  or  smaller.  Combining  the  veloc¬ 
ity  fields  to  produce  the  compromise  model  V5  allowed 
us  to  reduce  the  small-scale  noise  and  exploit  that  skill. 
This  result  is  consistent  with  the  results  obtained  by 
Murphy  (1988).  Based  on  both  numerical  results  and 
simple  theoretical  estimates,  Murphy  (1988)  showed 
that  an  “ensemble-mean  forecast  is  more  skilful  than  an 
individual  forecast”.  Spatial  smoothing  should  also  im¬ 
prove  the  skill  of  both  an  individual  and  an  ensemble- 
mean  forecast,  but  the  ensemble  mean  forecast  remains 
superior.  The  explanation  is  that  an  ensemble-mean 
forecast  provides  a  better  estimate  of  the  true  ocean 
circulation  than  a  single  one.  However,  because  of 
the  statistical  uncertainty,  these  conclusions  are  not 
definitive. 

From  a  practical  point  of  view,  one  would  like  to 
choose  the  best  model  to  make  forecasts.  Compar¬ 
ing  the  80th  percentile  and  RMS  trajectory  errors  in 
Figs.  8  and  9,  it  appears  that  V5  with  no  smoothing 
and  GLORYS  are  the  two  best  and  the  differences 


Compare  V5  and  optimally  smoothed  GLORYS 


Fig.  10  CDF  of  trajectory  error  as  in  Fig.  4  but  comparing  the 
optimally  smoothed  GLORYS  model,  red  line  (smoothed  with 
Gaussian  radius  of  about  167  km)  and  V5  model,  black  line 
(not  smoothed).  The  comparison  is  made  with  exactly  the  same 
set  of  drifter  trajectories.  Because  the  compromise  V5  model 
is  consistently  better  than  the  best  of  the  OGCMs  for  1-,  3- 
and  7-day  trajectory  estimates,  we  have  grounds  to  reject  the 
null  hypothesis  that  forming  the  compromise  model  improves 
estimates  only  because  it  has  smoother  velocity  fields 
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bclvvccn  them  arc  not  statistically  significant  at  the  95% 
confidence.  The  80lh  percentile  and  RMS  represent 
only  two  statistics  derived  from  the  full  distribution.  It 
is  informative  to  compare  the  full  distribution  by  com¬ 
paring  the  CDFs.  Figure  10  shows  the  CDF  of  trajectory 
error  for  GLORYS  optimally  smoothed  (R  =  160  km) 
and  V5  with  no  smoothing.  V5  is  consistently  better 
but  because  the  improvement  is  only  slight,  it  is  not 
statistically  significant  at  the  95%  level.  This  provides 
evidence  against  the  null  hypothesis  that  the  ensemble 
was  better  only  because  it  was  smoother.  Unfortunately 
with  the  limited  data  available,  we  cannot  reject  the 
null  hypothesis  with  95%  confidence.  One  factor  to 
bear  in  mind  is  that  the  weights  used  in  the  optimal 
model  combination  done  using  Eqs.  1  and  2  could  be 
a  function  of  the  filter  radius  used  in  the  smoothing, 
yet  the  weights  for  V5  where  chosen  from  the  RMS 
error  with  no  smoothing.  The  change  of  the  individual 
model  skill  with  the  spatial  smoothing  rather  supports 
this  hypothesis. 


6  Summary  and  discussion 

Ocean  state  estimation  with  data  assimilation  in 
OGCMs  with  sufficient  resolution  to  produce  an  ener¬ 
getic  mesoscalc  eddy  field  that  dominates  the  advection 
of  surface  drifters  is  a  very  young  science.  Adding  qual¬ 
ity  to  the  observational  data  is  clearly  a  great  challenge. 
Evidence  of  the  great  challenge  was  presented  herein, 
where  the  data  product,  SURCOUF,  revealed  more 
skill  than  the  five  individual  data-assimilativc  OGCMs 
in  estimating  the  trajectories  of  surface  drifters.  The 
OGCMs  all  assimilated  the  altimeter  data  upon  which 
the  SURCOUF  data  product  was  based.  Furthermore, 
our  results  were  for  the  equatorial  Atlantic  where 
SURCOUF  faced  the  additional  challenge  that  the 
geoslrophic  balance  breaks  down  within  about  5°  of  the 
equator.  In  this  version  of  SURCOUF,  the  calibration 
of  the  velocity  anomalies  was  completely  independent 
of  the  AOML  drifters  from  our  study  period.  However, 
recall  the  mean  currents  in  SURCOUF  were  dependent 
on  the  AOML  drifters  used  herein,  but  the  dependence 
was  barely  significant  at  the  95%  level,  see  discussion  in 
the  “Appendix”.  Another  important  fact  is  that  SUR¬ 
COUF,  as  a  data  product,  is  inherently  limited  to  esti¬ 
mates  of  past  events,  while  in  principle  the  OGCMs  run 
in  an  operational  configuration  can  make  predictions  of 
the  future,  although  we  did  not  assess  future  predictions 
herein. 

We  found  median  1-day  trajectory  errors  between 
about  15  and  20  km.  This  is  consistent  with  re¬ 
sults  found  with  the  US  Navy  EAS16  model  in  the 


East  China  Sea  (Huntley  ct  al.  2011).  Barron  et  al. 
(2007)  found  that  trajectory  estimation  accuracy  was 
regionally  dependent.  Although  their  NCOM  Equato¬ 
rial  Atlantic  results  (16  to  18  km)  fall  within  the  range 
we  found  for  all  models,  they  found  better  results  with 
NCOM  than  we  did,  presumably  because  their  region 
extended  to  ±15°  and  thus  included  more  quiescent 
regions.  Trajectory  error  of  course  increased  with  time 
and  by  7  days  we  found  median  errors  between  about  80 
and  115  km,  also  broadly  consistent  with  Barron  ct  al. 
(2007). 

In  practice,  controlled  experiments  with  eddying 
OGCMs  with  data  assimilation  are  extremely  costly 
in  terms  of  both  researcher  time  and  computational 
resources  and  need  to  be  well  justified.  Much  can  be 
learned  from  the  inter-comparison  of  existing  runs,  and 
the  results  described  herein  motivate  further  study. 
Our  results  of  accuracy  of  estimated  drifter  trajectories 
begs  the  question  why  SURCOUF  had  more  skill  than 
the  OGCMs  especially  since  all  the  OGCMs  assimi¬ 
late  the  altimeter  data  used  to  construct  SURCOUF. 
While  the  two  NRL  OGCMs  both  used  Navy  Oper¬ 
ational  Global  Atmospheric  Prediction  System  (NO¬ 
GAPS)  winds,  the  three  Mercator-Ocean  OGCMs  all 
used  ECMWF  winds,  as  did  SURCOUF.  This  raises 
questions  as  to  whether  the  temporal  resolution  of  the 
winds  is  an  important  factor  since  SURCOUF  used 
the  six-hourly  Interim  winds  from  ECMWF  while  the 
Mercator-Ocdan  OGCMs  were  forced  with  daily  aver¬ 
aged  operational  ECMWF  winds. 

We  speculate  that  some  benefit  could  be  obtained 
by  forcing  the  OGCMs  with  the  six-hourly  winds  rather 
than  daily  averages.  Indeed,  inertial  ocean  motions  are 
forced  by  moving  storms  when  the  wind  vector  rotates 
at  a  frequency  close  to  the  local  inertial  frequency  / 
(e.g.  Price  1981 , 1983).  The  inertial  period  27i/ /  is  about 
2.9  days  near  10°  N,  and  therefore,  the  rotation  of  the 
wind  stress  is  poorly  described  with  a  daily  forcing. 
Finally,  we  speculate  that  using  a  six-hourly  solar  forc¬ 
ing  may  have  stronger  consequences  than  using  a  six- 
hourly  wind  forcing  (as  compared  to  daily  forcings),  as 
found  by  B.  Barnier  (personal  communication)  in  the 
Mediterranean  Sea. 

None  of  the  models  included  tidal  forcing,  which  in 
the  deep  ocean  leads  to  errors  similar  to  that  due  to 
slippage  in  strong  winds,  about  1  cm/s.  However,  this 
is  unlikely  to  be  a  dominant  factor,  contributing  only 
about  a  kilometre  per  day  to  the  error. 

It  is  instructive  and  perhaps  surprising  that  one  of 
the  lowest  resolution  OGCMs,  GLORYS,  had  the  best 
skill  of  all  the  OGCMs.  The  fact  that  the  higher- 
resolution  PSY2  was  marginally  better  than  its  lower- 
resolution  counterpart  PSY3  suggests  that  resolution  is 
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by  no  means  a  hinderancc  and  perhaps  improves  the 
estimates  marginally.  Recall  from  Tables  2  and  3  that 
GLORYS  had  two  clear  advantages  over  its  counter¬ 
parts  at  Mercator-Ocdan:  The  Jason-1  and  ENVISAT 
altimeter  data  and  ARGO  and  other  T&S  data  were 
used  in  delayed  mode  instead  of  near  real  time.  De¬ 
layed  mode  allows  for  quality  control  with  removal  of 
outliers  and  other  suspect  data.  Furthermore,  GLO¬ 
RYS  integrated  the  incremental  analysis  update  (IAU) 
scheme  in  the  data  assimilation  method  (Bloom  ct  al. 
1996;  Bcnkiran  and  Greiner  2008),  which  is  believed 
to  be  an  improvement.  The  results  described  herein 
strongly  motivate  controlled  experiments  to  decisively 
determine  the  relative  importance  of  data  quality  and 
data  assimilation  scheme  in  producing  the  enhanced 
skill  of  GLORYS  over  PSY3  and  PSY2. 

Evidence  that  the  OGCMs  can  add  value  to  the 
data  was  obtained  from  the  result  that  compromise 
estimates  employing  all  the  OGCMs  had  skill  as  good 
or  better  than  SURCOUF.  Two  strategics  for  combin¬ 
ing  the  OGCM  results  to  produce  compromise  esti¬ 
mates  were  explored.  We  tried  averaging  the  latitude 
and  longitude  of  the  estimated  trajectory  from  each 
model,  weighted  by  its  mean  squared  error.  And  we 
tried  averaging  the  OGCM  velocity  fields  to  produce  a 
single,  compromise  velocity  field  that  was  used  to  make 
trajectory  estimates.  The  results  were  very  similar  for 
both  ensemble  methods  and  always  better  than  the  best 
single  model,  GLORYS;  compare  the  “ALL”  column 
of  Table  6  and  “GLORYS”  column  of  Table  4.  The  en¬ 
semble  trajectory  estimates  were  almost  always  better 
than  the  SLIRCOUF  estimates,  though  not  always  sta¬ 
tistically  significantly  so;  compare  the  “ALL”  column  of 
Table  6  and  “SURCF”  column  of  Table  4.  The  fact  that 
a  multi-model  ensemble  better  fits  the  observational 
signal  and  improves  the  skill  is  consistent  with  the 
results  from  Murphy  (1988)  where  the  superiority  of 
an  ensemble-mean  forecast  over  an  individual  forecast 
was  shown  both  with  numerical  ensemble  forecasts  and 
simple  theoretical  considerations.  The  result  from  the 
present  study  points  to  interesting  prospects  for  the 
users  of  operational  oceanography  forecasts. 

We  also  noted  that  most  OGCMs  had  an  “optimal 
low-pass  filter  radius”,  i.c.  removing  scales  smaller  than 
roughly  150  km  or  so  improved  the  estimates.  The  in¬ 
terpretation  is  that  the  OGCM's  surface  velocity  fields 
contain  errors  that  increase  in  importance  with  smaller 
scales.  Below  about  150  km,  the  errors  arc  more  im¬ 
portant  than  the  signal  (importance  assessed  by  impact 
on  Lagrangian  particle  trajectories).  These  errors  can 
be  attributed  to  various  sources.  Models  suffer  from 
shortcomings  like  their  lack  of  resolution,  errors  in  the 


forcing  fields  or  bulk  parametcrisations  used  which  all 
contribute  to  the  forecast  errors.  And  at  the  scales  con¬ 
sidered  (spatial  scales  less  than  the  Rossby  radius  and 
time  scales  less  than  7  days),  the  nonlinearities  of  the 
advection  terms  in  the  momentum  equations  are  more 
important  than  at  larger  scales  and  are  able  to  strongly 
amplify  small  initial  errors.  Furthermore,  all  numerical 
models  include  data  assimilation  where  the  model  state 
is  corrected  sequentially.  The  analysis  may  introduce 
some  shocks  when  the  increments  arc  applied,  i.c.  the 
thermodynamics  may  be  slightly  unbalanced,  especially 
at  the  equator  where  geostrophic  balance  does  not 
hold.  This  may  introduce  some  noise  in  the  surface 
currents,  especially  when  incremental  analysis  update 
is  not  used  (Ourmieres  et  al.  2006).  Lastly,  the  obser¬ 
vational  network  is  certainly  not  sufficient  to  correctly 
constrain  the  full  dynamics.  All  these  factors  contribute 
to  small-scale  errors.  By  filtering  the  small  scales,  one 
removes  those  errors  and  one  obtains  an  ocean  state 
estimation  closer  to  the  true  ocean.  As  a  consequence, 
the  trajectory  forecast  skill  is  improved.  This  result  is 
also  in  agreement  with  Murphy  (1988)  who  showed  that 
spatial  filtering  improves  both  ensemble-mean  and  in¬ 
dividual  forecasts.  Our  result  is  also  partially  consistent 
with  the  results  of  Huntley  et  al.  (2011)  who  found 
that  coarsening  the  US  Navy  EAS16  model  velocity 
fields  in  space  (or  in  time)  by  up  to  a  factor  of  8  via 
subsampling  did  not  degrade  the  trajectory  estimates. 
However,  buried  in  those  errors,  there  appeared  to  be 
some  useful  signal.  This  was  suggested  by  the  fact  that 
the  optimal  filter  for  the  compromise  velocity  field  that 
averaged  the  five  models  as  in  Eq.  2  was  between  0  and 
about  75  km  and  thus  smaller  than  the  100-  to  200-km 
optimal  filter  of  the  individual  OGCMs.  The  ensemble 
velocity  fields  may  have  had  skill  on  smaller  scales. 
Furthermore,  ensemble  velocity  fields  had  greater  skill 
than  even  the  best  single  model  optimally  smoothed, 
see  the  CDF  of  trajectory  error  in  Fig.  10.  But  because 
the  improvement  is  only  slight,  it  is  not  statistically 
significant  at  the  95%  level. 
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Appendix:  Models 

Description  of  PSY2  model  runs 

Mercator-Ocdan  PSY2  is  an  operational  ocean  analy¬ 
sis  and  forecasting  system  based  on  the  NATL12 
ocean  configuration  and  a  reduced  order  Kalman 
filter  data  assimilation  scheme.  The  NATL12  ocean 
model  configuration  is  a  regional  implementation  of  the 
occan/sca-icc  NEMO  numerical  framework  (Madec 
2(X)8)  carried  out  by  Mcrcator-Oc<5an  and  the  Euro¬ 
pean  DRAKKAR  collaboration  (DRAKKAR  Group 
2007).  The  configuration  is  similar  in  many  points  to  the 
one  used  by  Barnier  et  al.  (2006)  except  that  the  resolu¬ 
tion  is  1/12°  (Mercator  grid)  and  that  the  geographical 
domain  is  limited  to  the  North  Atlantic  (20°  S-70°  N) 
and  the  Mediterranean  Sea.  The  model  has  50  layers 
in  the  vertical  (Z  coordinate),  with  a  layer  thickness  of 
1  m  near  the  surface,  increasing  with  depth  (10  m  at 
50  m  depth)  up  to  450  m  at  5,500  m  depth.  It  uses  a 
partial  step  representation  of  the  bottom  topography 
and  a  momentum  advcction  scheme  that  both  yielded 
significant  improvements  (Penduff  et  al.  2007).  Parame- 
terisations  include  a  Laplacian  mixing  of  temperature 
and  salinity  along  isopycnals,  a  horizontal  biharmonic 
viscosity  and  a  turbulence  closure  scheme  (TKE)  for 
vertical  mixing.  The  turbulent  surface  fluxes  arc  com¬ 
puted  using  the  CLIO  bulk  formulation  (Goossc  et  al. 
2001)  forced  by  ECMWF  operational  analyses  (1-day 
averages).  The  surface  wind  stress  is  the  one  provided 
by  ECMWF  analyses  (i.e.  it  is  not  computed  using  the 
CLIO  formulation).  The  ocean  model  includes  the  river 
runoff  climatology  of  Dai  and  Trenbcrth  (2002). 

The  data  assimilation  scheme  is  based  on  the  singular 
evolutive  extended  Kalman  filter  formulation  proposed 
by  Pham  et  al.  (1998).  This  approach  has  been  used  for 
several  years  at  Mercator-Ocdan  and  was  implemented 
in  different  ocean  (re)analysis  systems  like  PSY2,  PSY3 
(http://bulletin.mercator-occan.fr/html/wclcome_cn.jsp) 
or  GLORYS  (see  below).  Details  about  the  implement¬ 
ation  of  the  data  assimilation  scheme  are  described  by 
Tranchant  et  al.  (2008).  A  key  aspect  of  the  method 
is  the  use  of  a  large  number  of  model  anomalies  (a 
few  hundred)  to  model  explicitly  the  background 
model  error  covariance.  In  PSY2.  the  control  vector 
consists  of  the  barotropic  height  and  the  temperature 
and  salinity  fields.  In  order  to  produce  a  balanced, 
analysed  ocean  state  the  velocity  is  deduced  from  the 
barotropic  sea  surface  height  and  mass  field  increments 
using  appropriate  physical  balance  operators.  Among 
them,  a  linear  equation  of  baroclinic  motion  is  used 
to  build  the  baroclinic  velocity  increments  near  the 


equator  (sec  Bcnkiran  and  Greiner  (2008)).  The  length 
of  the  assimilation  cycle  is  7  days,  and  the  increment 
is  applied  directly  to  the  model  state  at  the  analysis 
time.  The  assimilated  data  is  sea  level  anomaly  (in 
conjunction  with  the  mean  dynamic  topography  of  Rio 
and  Hernandez  (2004)),  SST  and  in  situ  profiles  (see 
Table  2). 

Description  of  PS  Y3  model  runs 

Mercator-Oc<5an  PSY3  is  very  similar  to  PSY2  (see 
above)  except  for  the  following  principal  differences. 
The  spatial  resolution  was  coarser  (1/4°  Mercator 
grid).  Furthermore,  the  error  covariance  statistics  were 
different  since  they  were  computed  from  a  reference 
experiment  with  no  data  assimilation  for  each  model 
configuration.  Finally,  the  SST  observation  operator 
was  slightly  different,  with  less  smoothing  of  the  model 
equivalent  for  PSY3,  and  slightly  different  weights 
applied. 

Description  of  the  GLORYS  model  runs 

GLORYS  is  a  project  with  objective  to  produce  a 
scries  of  realistic  (i.e.  close  to  the  existing  observations 
and  consistent  with  the  physical  ocean)  eddy  resolving 
global  ocean  rcanalyses.  The  version  1  of  stream  1 
(called  GLORYS1V1)  covering  the  Argo  years  (2002- 
2008)  is  used  in  this  study.  The  ocean  reanalysis  system 
used  in  GLORYS  is  similar  in  many  points  to  PSY2. 
However,  some  important  differences  exist  that  we 
review  in  the  following.  The  OGCM  used  in  GLO- 
RYS1V1  is  also  based  on  the  occan/sca-icc  NEMO  nu¬ 
merical  framework  (Madec  2008).  The  main  differences 
with  PSY2  arc:  (a)  the  configuration  is  global  (— 77S 
to  the  North  Pole)  and  (b)  the  horizontal  resolution  is 
coarser  (1/4°  Mercator  grid).  GLORYS1V1  shares  the 
same  vertical  grid  as  PSY2.  The  physical  parameterisa- 
tions  used  arc  basically  the  same,  except  that  GLORYS 
has  an  additional  harmonic  diffusion  operator  in  the 
tropical  band  (5.5°  N/5.5°  S)  to  improve  the  tropical 
instability  wave  physics  at  1/4°  horizontal  resolution. 
The  data  assimilation  scheme  is  close  to  PSY2,  except 
that  the  control  vector  includes  two  more  variables: 
the  zonal  and  meridional  velocity  components.  This 
allows  producing  a  velocity  increment  consistent  with 
the  analysed  mass  field  and  model  physics.  Because 
there  is  no  simple  relationship  between  the  mass  field 
and  the  circulation  near  the  equator  (e.g.  Bcnkiran  and 
Greiner  2008),  the  analysed  velocity  near  the  equator 
is  only  partially  applied.  The  velocity  increments  are 
set  to  zero  at  the  equator  and  increase  smoothly  with 
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latitude  to  become  maximal  at  7°.  The  last  important 
difference  with  PSY2  is  the  use  of  an  IAU  initialisation 
procedure  (Bloom  ct  al.  1996)  to  apply  the  increment 
to  produce  a  lime  continuous  ocean  analysis. 

Description  of  NCOM  model  runs 

The  2007-2008  NCOM  results  are  from  the  global 
ocean  forecast  system  (GOFS)  version  2.6  running  op¬ 
erationally  at  the  Naval  Oceanographic  Office.  Global 
NCOM  (Barron  ct  al.  2006)  uses  a  rotated  bipolar 
curvilinear  grid  with  spacing  about  20  km  at  the  equa¬ 
tor,  14  km  (1/8°)  at  midlatitudcs  and  5  km  in  the  Arctic. 
It  employs  an  implicit-free  surface  and  hybrid  vertical 
grid  with  19  terrain-following  sigma  levels  above  21 
fixed-level  z  coordinates  below  141  m,  chosen  to  be 
near  the  depth  of  a  typical  shelf  break.  The  layers  arc 
logarithmically  stretched  below  a  1-m-thick  uppermost 
layer.  In  GOFS  2.6,  data  assimilation  (Barron  ct  al. 
2007)  relaxes  to  a  background  of  Modular  Ocean  Data 
Assimilation  System  (Fox  et  al.  2002)  synthetic  profiles 
of  temperature  and  salinity  derived  based  on  remote 
observations  of  sea  surface  height  and  sea  surface  tem¬ 
perature  and  modified  by  in  situ  observations  through 
the  Navy  Coupled  Ocean  Data  Assimilation  System 
(NCODA;  (Cummings  2005)).  GOFS  obtains  wind 
stress  and  surface  heat  flux  forcing  from  the  NOGAPS 
(Rosmond  ct  al.  2002).  River  inflow  is  derived  from  a 
monthly  climatology  for  981  of  the  world’s  largest  rivers 
(Barron  and  Smedstad  2002).  'Tides  are  not  included  in 
the  present  global  model  and  are  optionally  added  in  a 
separate  step  when  nesting  a  higher-resolution  regional 
model.  Simulated  drifter  evaluations  with  earlier  ver¬ 
sions  of  global  NCOM  were  performed  by  Barron  et  al. 
(2007)  and  van  Scbille  et  al.  (2009). 

Description  of  HYCOM  model  runs 

IIYCOM  is  widely  used  by  the  ocean  community 
(http://www.hycom.org)  and  is  the  backbone  of  the 
global  eddy-resolving  (1/12°  horizontal  resolution)  real¬ 
time  nowcast/forccast  system  at  the  Naval  Oceano¬ 
graphic  Office.  The  model  has  nominal  1/12°  Mercator 
grid  horizontal  resolution  (6.5-km  grid  at  midlatitudes 
with  a  bipolar  patch  north  of  47°  N,  i.e.,  3.5-km  grid 
spacing  at  the  North  Pole)  and  32  hybrid  layers  in  the 
vertical  (pressure  coordinates  are  used  in  the  mixed 
layer,  isopycnal  coordinates  arc  used  in  the  ocean  in¬ 
terior  and  terrain-following  coordinates  are  used  in 
shallow  areas)  (Blcck  2002;  Chassignct  ct  al.  2003, 2006, 
2009). 

The  simulation  used  for  the  analysis  here  was 
restarted  from  a  spun-up  state  of  a  non-assimilative 


global  1/12°  HYCOM  simulation  with  climatological 
forcing.  The  assimilative  hindcast  began  in  May  2007 
and  integrated  through  the  analysis  period  while  be¬ 
ing  forced  by  the  three-hourly  NOGAPS  (http://www. 
nrlmry.navy.mil/nogaps_his.htm)  wind  stress,  wind 
speed,  heat  flux  (using  bulk  formula)  and  precipitation. 
Runoff  from  986  rivers  was  included  as  virtual  salinity 
flux  with  no  mass  exchange.  The  1/12°  assimilative 
global  HYCOM  system  has  been  validated  by  Metzger 
ct  al.  (2008, 2010). 

The  data  assimilation  in  HYCOM  is  performed  using 
the  NCODA  system.  NCODA  is  a  multivariate  optimal 
interpolation  scheme  that  assimilates  surface  obser¬ 
vations  from  satellites,  including  altimeter  and  multi¬ 
channel  sea  surface  temperature  data,  sea  ice  concen¬ 
tration  and  also  profile  data  such  as  XBTs,  CTDs  and 
profiling  floats  (Cummings  2005). 

Description  of  SURCOUF 

The  SURCOUF  currents  were  computed  globally  at 
a  six-hourly  temporal  resolution  using  the  following 
methodology  that  combines  estimates  of  gcostrophic 
current  and  the  Ekman  component  of  the  agcostrophic 
current.  Daily  gcostrophic  currents  were  obtained  from 
the  gradient  of  the  global  multi-mission  altimetric  maps 
distributed  by  AVISO.  The  mean  dynamic  topography 
(MDT)  used  to  reference  the  sea  level  anomalies  was 
the  recent  solution  computed  by  Rio  et  al.  (2011)  using 
the  following  three-step  method  (Rio  and  Hernandez 
2004).  First,  a  large-scale  mean  dynamic  topography 
was  obtained  from  the  CLS01  altimetric  mean  sea  sur¬ 
face  (Rio  and  Hernandez  2004)  and  a  geoid  model 
computed  at  Groupe  de  Researches  de  Geodesie  Spa- 
tiale,  Toulouse  from  4.5  years  of  GRACE  data.  Then, 
in  situ  measurements  (drifting  buoy  velocities  from 
1993  to  2008,  CTD  casts  and  ARGO  profiles  from 
1993  to  2007)  were  combined  with  altimetric  anomalies 
to  compute  synthetic  estimates  of  the  MDT  and  the 
corresponding  mean  currents.  The  synthetic  estimates 
were  finally  used  to  improve  the  large-scale  solution 
through  a  multivariate  objective  analysis. 

Six-hourly  ERA  INTERIM  wind  stress  fields 
(Berrisford  et  al.  2009)  were  used  to  compute  maps 
of  Ekman  currents  /7ek  through  the  use  of  a  simple 
statistical  model: 

wCk  =  exp  (/#)?, 

where  and  0  have  been  obtained  analysing  15-m 
drogued  drifting  buoy  velocities  that  were  collected  in 
the  framework  of  the  international  Global  Drifter  Pro¬ 
gram,  quality-controlled  and  distributed  by  the  AOML 
center.  To  estimate  //<*,  absolute  altimetric  velocities 
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Fig.  11  /I  (left)  and  9  (right) 
parameters  of  the  Ekman 
model  used  in  the 
computation  of  the 
SURCOUF  currents. 
Parameters  arc  displayed  by 
month  (from  /  Jan  to  12  Dec) 
and  latitude  (ranging  between 
10°  S  to  10°  N).  The  units  are 
(metres  per 

second)/(Newlons  per  square 
metre)  and  degrees, 
respectively 


1  2  3  4  5  6  7  8  9  1011  12 


-60-40-30-20-10  0  10  20  30  40  60 


were  interpolated  along  the  drifting  buoy  trajectories 
and  subtracted  from  the  buoy  velocities.  The  residual 
ageostrophic  current  was  further  filtered  using  a  1.25- 
to  20-day  band  pass  filter  to  focus  on  the  frequencies 
where  the  coherency  between  the  wind  stress  and  the 
Ekman  currents  is  maximal  (Rio  and  Hernandez  2003). 
Then  ERA  INTERIM  wind  stress  values  were  interpo¬ 
lated  along  the  drifting  buoy  trajectories  and  also  band¬ 
pass  filtered.  A  least  square  fit  was  finally  performed 
between  /}ek  and  f  so  as  to  obtain  the  and  0  parame¬ 
ters  by  latitudinal  bands  and  by  month.  The  resulting 
parameters  are  plotted  in  Fig.  1 1 .  For  this  specific  study, 
the  least  square  fit  was  performed  on  drifting  buoy 
velocities  available  for  the  year  2006  in  order  to  be 
independent  from  the  analysis  time  period  (2007-2008). 
Then  Ekman  currents  were  estimated  and  added  to 
the  geostrophic  currents  to  obtain  an  estimate  of  the 
total  surface  current.  Note  that  removing  the  drifters 
from  our  study  period  (July  2007  through  December 
2008)  from  the  SURCOUF  mean  fields  degraded  the 
trajectory  estimates  in  all  cases.  But  the  increase  in 
trajectory  error  was  similar  to  the  half- width  of  the  95% 
confidence  limits.  For  example,  the  80th  percentile  for 
3-day  trajectories  changed  from  70  ±  3  to  73  ±3  km 
in  Table  4,  while  for  7-day  trajectories  changed  from 
145  db  10  to  153  ±  10  km.  In  other  words,  the  change 
was  barely  significant  at  the  95%  level.  This  was  to  be 
expected  because  the  AOML  drifters  from  our  study 
period  and  main  analysis  region  were  not  the  dominant 
determinant  of  the  mean  field  there  since  less  than  15% 


of  the  244,608  drifter  velocity  measurements  within 
the  main  analysis  region  fell  within  the  July  2007  to 
December  2008  time  window. 
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